import numpy as np

def TVD_Limiter(uh, Nx, hx, dimPk, bcL, bcR, M=0.0):
    """实现TVD (Total Variation Diminishing) 限制器，用于数值计算中抑制激波捕捉时的非物理振荡
    
    该函数通过对每个网格点的多项式系数进行修正，确保数值解满足总变差不增性质，
    适用于有限体积法或间断伽辽金方法中的高阶格式处理。
    
    Args:
        uh (numpy.ndarray): 输入的多项式系数数组，形状为 (Nx, dimPk)
                            其中每行对应一个网格点的多项式系数，列对应不同阶数的系数
        Nx (int): x方向的网格单元数量，决定了需要处理的空间离散点总数
        hx (float): x方向的网格间距（单元长度）
        dimPk (int): 分片多项式空间P_k的维数，与多项式次数k相关（如k=2时dimPk=3），
                     决定了每个网格点的系数数量
        bcL (int): 左边界条件标识（如1表示周期边界，2表示常数值边界）
        bcR (int): 右边界条件标识（含义同bcL，用于右侧边界处理）
        M (float): 限制器强度参数，控制限制器的效果，M→∞相当于没有加限制器，M->0相当于TVD限制器
    
    Returns:
        numpy.ndarray: 经过TVD限制器修正后的多项式系数数组，形状与输入uh一致，
                       确保数值解在保持精度的同时抑制非物理振荡
    """
    uhb = np.vstack([[0, 0, 0], uh, [0, 0, 0]])
    
    # 初始化修改后的数组
    uhmod = np.zeros((Nx, dimPk))
    uhmod[:, 0] = uh[:, 0]  # 第一列保持不变
    
    # 设置边界条件
    if bcL == 1:
        uhb[0, :] = uh[-1, :]    # 左幽灵 = 最后一个物理单元
    elif bcL == 2:
        uhb[0, :] = uh[1, :]     # 左幽灵 = 第一个物理单元
    if bcR == 1:
        uhb[-1, :] = uh[0, :]    # 右幽灵 = 第一个物理单元
    elif bcR == 2:
        uhb[-1, :] = uh[-2, :]   # 右幽灵 = 最后一个物理单元
    
    # 遍历每个网格点计算限制器
    for i in range(Nx):
        # 计算deltaUR和deltaUL（注意MATLAB列索引1-based转Python 0-based）
        deltaUR = uh[i, 1] + (2/3) * uh[i, 2]
        deltaUL = uh[i, 1] - (2/3) * uh[i, 2]
        
        # 计算相邻网格的差值
        deltaURM = uhb[i + 2, 0] - uhb[i + 1, 0]
        deltaULM = uhb[i + 1, 0] - uhb[i, 0]
        
        # 应用minmod限制器
        deltaURM1 = minmod(deltaUR, deltaURM, deltaULM, hx, M)
        deltaULM1 = minmod(deltaUL, deltaURM, deltaULM, hx, M)
        
        # 更新修改后的数组
        uhmod[i, 1] = (deltaURM1 + deltaULM1) / 2
        uhmod[i, 2] = 3 * (deltaURM1 - deltaULM1) / 4
    
    return uhmod


def minmod(a, b, c, hx, M=0.0):
    """实现minmod（最小模）限制器核心逻辑，用于TVD格式中选择合理的梯度/差值，抑制非物理振荡
    
    该函数通过判断输入差值的符号一致性和数值大小，选择“最小模”的有效差值，
    避免因梯度过度估计导致的数值解振荡
    
    Args:
        a (float): 待判断的目标差值/梯度（如当前网格点的多项式系数推导的梯度）
        b (float): 参考差值1（通常为当前网格与右侧相邻网格的物理量差值）
        c (float): 参考差值2（通常为当前网格与左侧相邻网格的物理量差值）
        hx (float): x方向的网格间距（网格单元的长度），用于计算数值阈值，避免小量干扰
        M (float): M → ∞相当于没有加限制器. M->0 相当于TVD限制器
    
    Returns:
        float: 经过minmod规则筛选后的有效差值/梯度，规则如下：
               1. 若a的绝对值小于M*hx²（小值阈值），直接返回a（避免过度修正）；
               2. 若a、b、c符号一致，返回与它们同符号且绝对值最小的值（选择最小模梯度）；
               3. 若a、b、c符号不一致，返回0（抑制可能引发振荡的异常梯度）
    """
    
    if abs(a) < M * (hx ** 2):
        a1 = a
    else:
        # 检查三个值的符号是否一致
        sign_a = np.sign(a)
        sign_b = np.sign(b)
        sign_c = np.sign(c)
        
        if sign_a == sign_b and sign_a == sign_c:
            # 取绝对值最小值并保持符号
            a1 = sign_a * min(abs(a), abs(b), abs(c))
        else:
            a1 = 0
    
    return a1
    